rm(list=ls())
library(readstata13)
library(plyr)

setwd("")

confidence_interval <- function(vector, interval) {
  vec_sd <- sd(vector, na.rm = T)
  n <- length(vector[!is.na(vector)])
  vec_mean <- mean(vector, na.rm = T)
  error <- qt((interval + 1)/2, df = n - 1) * vec_sd / sqrt(n)
  result <- c("lower" = vec_mean - error, "upper" = vec_mean + error)
  return(result)
}

d <- read.dta13("prediction_fig1.dta")

after <- NULL
for (i in 1:length(unique(d$id))){
  after_coup <- d$xb[d$id==unique(d$id)[i] & d$post_fail==1]
  if (length(after_coup) > 0) {
    after[[i]] <- after_coup
  }
}
after <- after[lapply(after,length)>0]
lapply(after, length)

preds_after <- ldply(after, rbind)
t_preds_after <- apply(preds_after, 2, mean, na.rm=T)
ci_after <- t(apply(preds_after, 2, function(x) confidence_interval(x,.95)))

before <- NULL
for (i in 1:length(unique(d$id))){
  before_coup <- rev(d$xb[d$id==unique(d$id)[i] & d$post_fail==0]) #reverse 1 so we can add NAs to proper side
  if (length(before_coup) > 0) {
    before[[i]] <- before_coup
  }
}
before <- before[lapply(before,length)>0]
lapply(before, length)

preds_before <- ldply(before, rbind)
preds_before <- rev(preds_before) # Rev again so right probs are AT END of period, not start
t_preds_before <- apply(preds_before, 2, mean, na.rm=T)
ci_before <- t(apply(preds_before, 2, function(x) confidence_interval(x,.95)))

length_plot_x_each_side <- 100
cut_after <- t_preds_after[1:length_plot_x_each_side]
cut_before <- t_preds_before[(length(t_preds_before)-length_plot_x_each_side-1):length(t_preds_before)]
ci_95_upper_before <- ci_before[(length(t_preds_before)-length_plot_x_each_side-1):length(t_preds_before),2]
ci_95_lower_before <- ci_before[(length(t_preds_before)-length_plot_x_each_side-1):length(t_preds_before),1]
ci_95_upper_after <- ci_after[1:length_plot_x_each_side,2]
ci_95_lower_after <- ci_after[1:length_plot_x_each_side,1]
all <- c(cut_before, cut_after)
ci_upper <- c(ci_95_upper_before, ci_95_upper_after)
ci_lower <- c(ci_95_lower_before, ci_95_lower_after)

#PLOT
years = 5.5
months = years*12/2
range = (102-months):(102+months)

#pdf("/Analysis/Figures/Figure1.pdf", width=20, height=13)
par(mar=c(7.5,5,3,2))
plot(NULL, NULL, xlim=c(0,12*years), ylim=c(min(all[range]), max(all[range])), 
     type="l", bty="l", ylab="", xaxt="n", cex.lab=3.5, cex.axis=3.5, xlab="")
lines(lowess(all[range], f=0.06), lwd=4, col="black")
abline(v=months, col="maroon", lwd=3)
axis(1, at=c(0, 12*years/6, 12*years/6*2, 12*years/6*3, 12*years/6*4, 12*years/6*5, 12*years), 
     labels = c(36, 24, 12, 0, 12, 24, 36), cex.axis=3.5, padj = 0.5)
mtext("Number of Months Before and After a Failed Coup", 1, line=6, cex=3.5)
#dev.off()

## GENERAL:
par(mar=c(5,5,3,2))
plot(NULL, NULL, xlim=c(0,length(all)), ylim=c(min(all), max(all)), 
     type="l", bty="l", ylab="", xaxt="n", cex.lab=2.2, cex.axis=2.2,
     xlab="Number of Months Before and After a Failed Coup", lwd=2)
lines(lowess(all, f=0.03), lwd=4, col="black")
abline(v=length(all)/2-1, col="maroon", lwd=3)
axis(1, at=c(0, 12*years/4, 12*years/4*2, 12*years/4*3, 12*years), 
     labels = c(24, 12, 0, 12, 24), cex.axis=2.2)





#### FIGURE 1:
library(readxl)
margins <- read_excel("Data/results_colpus.xlsx", col_names = c("a", "b", "c"))
margins <- data.frame(t(margins[c(1,5,6),]))
colnames(margins) <- margins[1,]
margins <- margins[2:3,]
margins <- apply(margins,2,  as.numeric)
pdf("Analysis/Figures/FigureA3_mareff.pdf", width=12, height=8)
par(mar=c(5,5,4,2))
plot(c(0, 1),margins[,1], xlim=c(-0.25, 1.25), ylim=c(min(as.numeric(margins[,2])), max(as.numeric(margins[,3]))), xaxt="n",
     ylab="", xlab="Failed Coup", lwd=2, col="blue", bty="l", pch=19, cex=2, , cex.axis=2,
     cex.lab=2,  main="COLPUS", cex.main=3)
axis(1,at=c(0, 1), cex.axis=2)
segments(0, margins[1,2], 0, margins[1,3], lwd=3, col="blue")
segments(1, margins[2,2], 1, margins[2,3], lwd=3, col="blue")
dev.off()


